# R script for figures in "Human Rights versus National Interests: Shifting US Public Attitudes on the International Criminal Court"
# Author: Kelebogile Zvobgo

# ORIGINAL CODE FROM:
# Chaudoin, Stephen, 2013, "Replication data for: Promises or Policies? 
# An Experimental Analysis of International Agreements and Audience Reactions", 
# https://hdl.handle.net/1902.1/21980, Harvard Dataverse, V2, UNF:5:0/Z3wzEVAt4tzNUW3vdp0w==; 
# replication_2013_07_29.R [fileName]



setwd("/Users/kelebogilezvobgo/Desktop/Research/1_ICC-Public-Opinion/ISQ_FINAL/Data")

library(foreign)
library(MASS)
library(gplots)
library(ggplot2)
library(Hmisc)
library(sciplot)
library(RItools)

rm(list=ls())

qtdata <- read.dta("HR-versus-NI_ISQ-Replication_Data.dta")
qtdata2 <- qtdata[ which(qtdata$treat_group!="NA" & qtdata$support_join_court !="NA"), ]
attach(qtdata2)

nbyti<- as.matrix(aggregate(support_join_court ~ treat_group, data=qtdata2, FUN=function(x) sum(!is.na(x)) ))

n0<-nbyti[1,2]
n1<-nbyti[2,2]
n2<-nbyti[3,2]
n3<-nbyti[4,2]

qtdata2_app <- qtdata2[ which(qtdata2$support_join_court =="1"), ]

abyti<- as.matrix(aggregate(support_join_court ~ treat_group, data=qtdata2_app, FUN=function(x) sum(!is.na(x)) ))

a0<-abyti[1,2]
a1<-abyti[2,2]
a2<-abyti[3,2]
a3<-abyti[4,2]

n<-5000

pi_0 <- as.matrix(rbeta(n, a0+0.5, n0-a0+0.5))
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
pi_3 <- as.matrix(rbeta(n, a3+0.5, n3-a3+0.5))

t0 <- as.matrix(rep(0,n, nrow=n, ncol=1))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))

tis <- as.matrix(rbind(t0,t1,t2,t3))
pis <- as.matrix(rbind(pi_0,pi_1,pi_2,pi_3))

betas <- as.data.frame(cbind(tis,pis))
betas$treat_group = factor(betas$V1, labels = c("Control", "Human Rights","National Interest","Competitive"))

###
# Figure ##
###
#This is the main plot, "Treatment Effects, All Respondents, DV = Join Court"

lineplot.CI(x.factor = betas$treat_group, response = betas$V2, type="p",
            data = betas, main = "", xlab = "Treatment Group", ylab = "Approval %",
            ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )